module retired_valuer_mod
USE functions_mod
USE parameters_mod
IMPLICIT NONE
CONTAINS


subroutine retired_valuer(capital,kprime_l,kprime_h,ss,tr,r,psi,vprime,v_out,k_out,avg_earnings,age,lambdak_in)
DOUBLE PRECISION,    intent(in)   :: capital,ss,tr,r,psi,avg_earnings,lambdak_in
DOUBLE PRECISION,    intent(in)   :: vprime(kgridN),kprime_l,kprime_h
INTEGER, intent(in)   ::age
DOUBLE PRECISION,    intent(out)  :: v_out,k_out
DOUBLE PRECISION                  :: v_out_temp,k_out_temp

v_out=-brent(kprime_l,capital,kprime_h,obj_utils,.00000001D0,k_out,capital,ss,tr,r,psi,vprime,avg_earnings,age,lambdak_in)

v_out_temp=-brent(kprime_l,max(min(k_out+.02D0,kprime_h-.001D0),kprime_l+.00001D0),kprime_h,obj_utils,.00000001D0,k_out_temp,capital,ss,tr,r,psi,vprime,avg_earnings,age,lambdak_in)
if(v_out_temp>v_out) then
v_out=v_out_temp
k_out=k_out_temp
end if

v_out_temp=-brent(kprime_l,max(min(k_out-.02D0,kprime_h-.001D0),kprime_l+.00001D0),kprime_h,obj_utils,.00000001D0,k_out_temp,capital,ss,tr,r,psi,vprime,avg_earnings,age,lambdak_in)
if(v_out_temp>v_out) then
v_out=v_out_temp
k_out=k_out_temp
end if


end subroutine retired_valuer


DOUBLE PRECISION FUNCTION obj_utils(kprime,capital,ss,tr,r,psi,vprime,avg_earnings,age,lambdak_in)
DOUBLE PRECISION,    intent(in)   :: capital,kprime,ss,tr,r,psi,lambdak_in
DOUBLE PRECISION,    intent(in)   :: vprime(kgridN),avg_earnings
INTEGER, intent(in)   ::age
DOUBLE PRECISION                  :: vprime_out(1),kprime_pass(1),con_tot,dummy,con,eng
    con_tot=retcon(capital,kprime,tr,r,ss,avg_earnings,lambdak_in)
    kprime_pass(1)=kprime
    CALL spline_linear_val1(kgridN,kgrid,vprime,kprime_pass(1),vprime_out(1),1)
    dummy=energy_sharer(0.00D0,con_tot*.15D0,con_tot,energy_share_equation,.0001D0,con,con_tot,age)
    eng=(con_tot-con)/(pe(year_calculating)+taue)
    obj_utils=-valuer_ret(con,eng,psi,vprime_out(1),adult_equivalence(age))

END FUNCTION obj_utils








DOUBLE PRECISION FUNCTION brent(ax,bx,cx,func,tol,xmin,capital,ss,tr,rate,psi,vprime,avg_earnings,age,lambdak_in)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
INTEGER, intent(in)   ::age
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol,capital,avg_earnings,lambdak_in
    DOUBLE PRECISION, INTENT(IN)    :: ss,tr,rate,psi,vprime(kgridN)
    DOUBLE PRECISION, INTENT(OUT)   :: xmin

INTEGER, PARAMETER          :: ITMAX=100
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,capital,ss,tr,rate,psi,vprime,avg_earnings,age,lambdak_in)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,capital,ss,tr,rate,psi,vprime,avg_earnings,age,lambdak_in)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent












subroutine spline_linear_val1 ( ndata, tdata, ydata, tval, yval, ndata2 )
      INTEGER :: left, right,i
      INTEGER, intent(in)  ::ndata,ndata2
      DOUBLE PRECISION,  intent(in) :: tdata(ndata),ydata(ndata),tval(ndata2)
      DOUBLE PRECISION,  intent(out) :: yval(ndata2)
      DOUBLE PRECISION  ::  ypval
ypval=0.0D0
yval=0.0D0
do i=1,ndata2
    call rvec_bracket1 ( ndata, tdata, tval(i), left, right )
  ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )
  yval(i) = ydata(left) +  ( tval(i) - tdata(left) ) * ypval
end do
  return
end SUBROUTINE spline_linear_val1

subroutine rvec_bracket1 ( n, x, xval, left, right )
  implicit none
      INTEGER, intent(in)  ::n
      INTEGER, intent(out)  ::left,right
      DOUBLE PRECISION, intent(in)   ::xval,x(n)
  integer i
  do i = 2, n - 1

    if ( xval < x(i) ) then
      left = i - 1
      right = i
      return
    end if
   end do
  left = n - 1
  right = n
  return
end subroutine rvec_bracket1


END module retired_valuer_mod


